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Abstract 

Energetic antiprotons in cosmic rays can serve as an important indirect signature of 
dark matter. Conventionally, the antiproton flux from dark matter decay or annihila- 
tion is calculated by solving the transport equation with a space-independent diffusion 
coefficient within the diffusion zone of the galaxy. Antiproton sources outside the dif- 
fusion zone are ignored under the assumption that they propagate freely and escape. 
In reality, it is far more likely that the diffusion coefficient increases smoothly with 
distance from the disk, and the outlying part of the dark matter halo ignored in the 
conventional approach can be significant, containing as much as 90% of the galactic 
dark matter by mass in some models. We extend the conventional approach to address 
these issues. We obtain analytic approximations and numerical solutions for antipro- 
ton flux for a diffusion coefficient that increases exponentially with the distance from 
the disk, thereby including contributions from dark matter annihilation/decay in es- 
sentially the full dark matter halo. We find that the antiproton flux predicted in this 
model deviates from the conventional calculation for the same dark matter parameters 
by up to about 25%. 



1 Introduction 



While astronomical observations have firmly established that about 80% of the matter con- 
tent of the universe exists in the form of non-baryonic dark matter, its microscopic prop- 
erties remain hitherto unknown. Weakly interacting massive particles (WIMPs) are the 
best-motivated candidates for dark matter from a theoretical point of view. WIMPs couple 
weakly to standard model (SM) particles, opening the possibility of dark matter particles 
decaying or annihilating into SM final states. Energetic cosmic rays produced in such decay 
or annihilation processes - in particular antimatter, which is rarely produced in astrophys- 
ical processes - can serve as important indirect signatures of dark matter in the galaxy. 
Observation of such signals can reveal information on the microscopic properties of dark 
matter. 

Our focus in this paper is on antiprotons produced in the annihilation or decay of dark 
matter in the Milky Way galaxy. This has become a topic of significant interest following 
recent measurements of the antiproton flux and the ratio of antiproton to proton flux up to 
180 GeV by the PAMELA experiment [1, 2}. No deviations from the expected astrophysical 
flux were observed. This data can in principle be used to put bounds on dark matter 
properties; however, to do this, the effects of propagation of antiprotons between the dark 
matter decay/annihilation location and the detector must be properly accounted for. 

Propagation of antiprotons in the galaxy is governed by transport equations, which in- 
corporate interactions with galactic magnetic fields and the interstellar medium (ISM). Con- 
ventional calculations solve the transport equation in a two-zone model: a thin disk of (9(100 
pc) thickness on the galactic plane where the interactions with the ISM occur, embedded in 
a larger diffusion region where galactic magnetic fields are appreciable and trap cosmic rays. 
The diffusion region is taken to be a cylinder of radius R of order 20 kpc and half-thickness 
L of order 1-10 kpc, and the diffusion coefficient is assumed to be position-independent in- 
side this region. Outside this diffusion region, antiprotons are assumed to propagate freely, 
leading to vanishing antiproton density on the boundaries of the cylinder in the steady state. 
With these assumptions, the transport equations can be solved in a straightforward manner. 
In particular, for antiproton energies of interest for dark matter searches (of order 10 GeV 
and above), additional approximations can be made that allow an analytic solution to the 
antiproton density and flux in the form of a Bessel series [3] . 

The conventional two-zone model is, however, a rather crude approximation. First, the 
assumption of a sharp boundary at L between the diffusion and free-propagation zones is 
unphysical. In reality, magnetic fields are known to decrease gradually with the distance from 
the disk, with exponential decay providing a reasonable fit to data. The diffusion coefficient 
likely follows a similar exponential profile [4] . Second, a typical dark matter halo is spherically 
symmetric and extends beyond the diffusion region and into the free propagation zone, in 
particular in the vertical direction: for example, for an isothermal dark matter profile in a 
model with L — 1 kpc, the diffusion zone contains only 10% of the dark matter mass of the 
full halo. The two- zone model completely ignores the antiprotons produced by dark matter 
decay/annihilation outside the diffusion zone. 
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The aim of this paper is to extend the conventional formalism for antiproton flux calcula- 
tions to overcome these shortcomings. 1 After reviewing the conventional formalism in Sec. 2, 
we consider a three- zone model in Sec. 3. This model still assumes an abrupt change in the 
diffusion coefficient at L, but includes the sources in the free-propagation zone extending 
to D > L, so that essentially all of the antiprotons from the dark matter halo are taken 
into account. We find an analytic solution to this model, and show that in the limit when 
diffusion in the free-propagation zone is completely absent, the sources in this zone have no 
effect on the flux measured at Earth. We then consider a model with an exponentially vary- 
ing diffusion coefficient [4] in Sec. 4, and obtain both a numerical solution and two analytic 
approximations which converge to it at high energies (above 50 GeV or so, depending on 
the desired accuracy). In Sec. 5, we present numerical results quantifying the effects of this 
more realistic propagation model on the antiproton fluxes produced by dark matter decay 
or annihilation, and the corresponding bounds on dark matter properties. We close with 
concluding remarks in Sec. 6. 

2 Antiproton Flux: The Conventional Formalism 

The starting point for conventional calculations is the steady-state diffusion equation for 
antiprotons [3]: 

- V [K(x, E) Vn p \ + ^-(V c (z)n p (E, x)) + 2h5(z)T ann n p (E, x) = q p (x, E) , (1) 

where n p is the antiproton density, K is the diffusion coefficient, the convective wind term 

V c (z) = sign (z)V c (2) 

is representative of motion of the medium responsible for diffusion, and q is the source term 
describing antiproton production in dark matter annihilation or decay. 2 The 2h8{z)Y ann term 
corresponds to antiproton interactions with the Interstellar Medium (ISM); h — 0.1 kpc is 
the half-width of the interaction zone, and T ann is the annihilation rate between antiprotons 
and protons, parameterized as [6] 

T ann = (n H + ^n He )o°™V p , (3) 

1 The extension presented here is similar to the formalism developed by us in a previous paper [5] for 
positrons. However, the absence of energy-loss terms for antiprotons allows us to make significant progress 
analytically, which was not possible for positrons. 

2 The source also contains the contribution of the so-called "tertiary" term, describing inelastic antiproton 
collisions with the ISM. We will ignore this contribution throughout this paper, as it is expected to be 
numerically unimportant for dark matter studies. In principle, it could be included in the calculations using 
an iterative approach. 
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Model 


6 


K (kpc 2 /Myr) 


L (kpc) 


V c (km/s) 


MIN 


0.85 


0.0016 


1 


13.5 


MED 


0.70 


0.0112 


4 


12 


MAX 


0.46 


0.0765 


15 


5 



Table 1: Galactic propagation models [9, 10]. All models set R = 20 kpc. 



ann/T-i\ 

a m U ) = 



i 



where nn ~ 1 cm 3 and n# e ~ 0.07n# are the hydrogen and helium number densities, and 
a|; n (T) is given by [7, 8] 

f 661(1 + 0.0115(T/GeV)-°- 774 - 0.948(T/GeV) a0151 ) mb, T < 15.5 GeV; 
36(T/GeV)~ a5 mb, otherwise. 

(4) 

Energy loss due to bremsstrahlung and collisions with cosmic microwave background photons, 
an important aspect of positron propagation, is a negligible effect for the more massive 
antiprotons and hence absent in the diffusion equation (1). Following [3, 15], we have also 
ignored solar modulation effects, which are expected to be unimportant at the high p energies 
(10 GeV and above) of interest for dark matter searches. 

The conventional approach is to assume an energy dependent but position independent 
diffusion coefficient 

K{E) = K p(n GY ) s , (5) 

where (3 — v/c and TZqv = pc/eZ is the rigidity of the particle measured in gigavolts (GV). 
Eq. (1) is then solved in a cylindrical region of radius R and half-thickness L, with vanishing 
antiproton density at the boundaries of the cylinder. The parameters L, R, K , Vc and 5 
define the galactic propagation model. We list three commonly used models in Table 1; these 
values are obtained from the analysis of observed isotope ratios in cosmic rays, primarily 
boron to carbon (B/C) ratio [10]. 

An analytic solution can be obtained by expanding Tip dbS db Bessel series 

n p (p,z,E) = Y,N i (z,E)J , (6) 

where Jo is the zeroth order Bessel function of the first kind, and Q's are the zeros of Jo- This 
reduces Eq. (1) to a set of ordinary differential equations on Ni(z), with E acting simply as 
a label. The source for each Ni is given by the Bessel transform of the source term, 

= f* dppJo ) q(p, z) , (7) 

where J\ is the first-order Bessel function of the first kind. The solution has the form [3] 
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Figure 1: Left: The dark matter halo extends significantly beyond the diffusion zone, but only 
sources inside the zone are considered in the conventional formalism (Section 2). Right: The 
three-zone formalism of Section 3 includes sources in the free propagation zone in addition 
to the diffusion zone. 

where we defined a = Vc/(2K), as well as 

*2 \ 1/2 



S, = 2\a 2 + C ' 



R 2 J KS % 

Vi{z) = 2 f Z e^- z 'hmh[S t (z-z')/2] qi (z')dz'. 
Jo 



A, = Vc + 2kTann - Bi = KS % [A, + coth^L/2)] 



(9) 



In particular, at the position of the Earth (corresponding to z = 0), where the fluxes are 
measured, the solution simplifies to 



e- aL yi {L) 



(10) 



(11) 



B iS mh(SiL/2) 

The antiproton flux at the top of the Earth's atmosphere can then be calculated as 

$_(£) = Eln p (r Q ,z Q = Q,E), 
where r = 8.5 kpc is the distance from the solar system to the galactic center. 

3 Three-Zone Propagation Model 

A straightforward way to include contributions from sources outside of the diffusion zone is to 
consider the transport equation in a cylinder of half-thickness D > L, choosing D sufficiently 
large so that all or virtually all of the dark matter halo is contained within this cylinder (see 
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Fig. 1, right). The diffusion coefficient remains of the form (5) for |^|<L, whereas in the free 
propagation region L < \z\ <Dwe assume 

K{E) = K^ilZGvY (12) 

with Ki ^> Kq, corresponding to significantly weaker magnetic fields and longer diffusion 
lengths in this region. To describe completely free propagation, one can take K\ — > oo 
at the end of the calculation. Note that the choice of the energy dependence in Eq. (12), 
which was taken to be the same as in the diffusion region, is a matter of mathematical 
convenience; we do not expect the results to depend on this choice in the large K\ limit. We 
impose the boundary condition n p = on the boundaries of the extended cylinder, since the 
(extragalactic) sources and magnetic fields outside the cylinder can be safely neglected. 

The transport equation in the three- zone approach can still be solved analytically. Ex- 
panding np in Bessel series as before, the solutions for Ni inside each zone are obtained 
as in the conventional approach. For each JVj, the solution in each zone contains two free 
constants parametrizing the solutions of the homogeneous (sourceless) equation. Two of 
these constants are determined by boundary conditions (iVj = at z = D and the matching 
condition at z — 0, which is the same as in the conventional approach). The other two are 
obtained from the matching conditions at the boundary between the zones at z = L: 

N t (L-e)=N l (L + e), K^{L - e) = K^{L + e) . (13) 

dz dz 

Physically, these conditions describe the continuity of antiproton density and flux, respec- 
tively, across the boundary; the latter condition is independent of the convection term, which 
is continuous at z — L. The resultant solution in the diffusion zone (\z\ < L) has the form 

Ni(z) = e a (l*l-^ [cosh(^/2) + A sinh(^/2)] - ^ . (14) 

Pi -ft &i 

The numerator and denominator in the coefficient are given by 

a t = UD)e- 2 ' aA + l^sinh^A) + — ^(L) cosh(^A) ; 

(5i = KSilAiSmhiSiL/2) + cosh(^L/2)] cosh(^A) 

+KSi [Ai cosh(SiL/2) + sinh(^L/2)] sinh(^A) . (15) 

Here, Si, Ai, a, and i/i are as defined in Section 2, and we also defined 

5 = ^, S i = 2\a' + ±) , A = (D-L)/2, 

Yi{z) = 2 T e a ^- z ^cosh{S l {z-z')/2)q l {z')dz', 
Jo 

Uz) = 2£e^ z - z 'hmh(S l (z-z')/2)q t (z')dz'. (16) 
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In particular, the flux at Earth is given by 

■^-^ 

From the above solution, it can be easily seen that in the limit K — > oo the solution 
inside the diffusion zone reduces to the one obtained in the conventional formalism, up to 
corrections of order K/K. In other words, the sources outside the diffusion zone do not 
contribute to the flux observed on Earth: all antiprotons from those sources get reflected by 
the boundary of the diffusion zone and do not penetrate it. This result in a sense justifies the 
use of conventional formalism, which ignores such sources. On the other hand, this result 
is only valid assuming an infinite, and infinitely sharply localized, jump in the diffusion 
coefficient at the boundary of the diffusion zone; both assumptions are clearly unphysical. It 
is far more reasonable to expect that the diffusion coefficient changes smoothly on the few- 
kpc length scale away from the disk. This motivates the analysis of the following section. 



4 Exponentially Increasing Diffusion Coefficient 

Since diffusion is caused by charged particles getting confined by galactic magnetic fields, 
variations in the diffusion coefficient should follow variations in the magnetic field strength. 
The magnetic field in the galaxy, while not precisely known, is believed to follow the approx- 
imate profile [11] 

^^PW^fe-^)' < 18 > 

The diffusion coefficient is expected to have a similar exponential spatial dependence. Ig- 
noring radial dependence, we can model the diffusion coefficient as [4]: 

K(E, z) = K e (E) exp(|z|M), K e (E) = K 2 p {^f-^ . (19) 

A numerical code for propagation of cosmic rays with a diffusion coefficient of this form 
has been presented in Ref. [4]. Using this code, it was demonstrated that a consistent 
fit to observed cosmic ray isotope ratios can be obtained in this model. Typical values 
of the parameters producing consistent fits are as follows: z t = 4 kpc, 5 = 0.57, K 2 = 
(0.55 x 10 28 cm 2 s _1 kpc _1 ) x z t = 2.2 x 10 28 cm 2 s~ 1 . We will use these values in all numerical 
calculations and plots in this paper. 

Our goal is to solve Eq. (1) with the diffusion coefficient of the form (19) and a boundary 
condition of vanishing antiproton density at large distances away from the galaxy. In practice, 
we demand rip = at \z\ = L and p = R, choose L and R sufficiently large so that 
essentially all of the dark matter halo is contained in the cylinder \z\ < L, p < R, and solve 
the equation within this cylinder. Since the setup still possesses cylindrical symmetry, the 
Bessel transform (6) can be used as before to reduce the transport equation to an (infinite) 
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set of ordinary differential equations. We were unable to find closed-form analytic solutions 
of these ODEs for the diffusion coefficient of the form (19) and the wind term of the form (2), 
but it is straightforward to solve them numerically, terminating the Bessel series after a finite 
number of terms, Nb- We verified that with typical dark matter profiles, the Bessel series 
converges rapidly, and convergence at a 1% level is achieved for Nb < 50 in all cases we 
studied (see Section 5). We also find it useful to consider two situations, slightly different 
from the real one, in which an analytic solution can be easily found: 

• Case 1: The convective wind term has the same exponential profile as the diffusion 
coefficient, namely V c (z) = V e exp(\z\/z t ); 

• Case 2: The convective wind term is neglected. 

While both situations are unphysical, they in a sense "bracket" the desired one (constant 
non-zero wind term): The first one underestimates the antiproton flux due to the artificially 
large wind term at large z, while the second one overestimates the flux. For case 1, the 
solution closely resembles Eq. (8): 

= e " M ~ L) BJmL/2) WW) + A, sinh( SiZ /2)] - fg , (20) 



although with slightly different definitions: 



v - 1 a 



2K e 2z t > \ R 2 z t K, 




1/2 



A = Ve + + "4; *i = K eSi [A + cothOSL/2)] , 

yAz) = 2 T e a(z ' z,) sinhlSAz - z')/2} qAz') e~ z ' ' Zt dz' . (21) 
Jo 

For case 2, one can simply set V e = in the above solution. As the wind term becomes 
negligible at high antiproton energies, the case 1 and 2 fluxes should approach each other, 
and therefore the true solution, in this regime. As we will see in the next section, the 
antiproton fluxes predicted in cases 1 and 2 are actually quite close to each other throughout 
the energy range relevant for dark matter searches. 



5 Results 

We model the dark matter distribution in the Milky Way galaxy with two widely used 
profiles, isothermal and Einasto. The isothermal profile is given by 

P(r) = Pe ^f'L (22) 
1 + [r/r s y 
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Figure 2: Ratio of antiproton flux computed from exponential diffusion coefficient model to 
flux computed from conventional model for isothermal (top) and Einasto (bottom) profiles. 
The left (right) panels show the results for decaying (annihilating) dark matter. Dashed 
and dot-dashed lines correspond to fluxes analytically calculated with and without the ex- 
ponential wind term respectively. The solid line in between corresponds to fluxes calculated 
numerically with a constant convective wind term. 



where p & = 0.3 GeV cm 3 is the local dark matter density in the solar neighborhood, and 
r s = 5 kpc. The Einasto profile is [12] 



p(r) = p& exp 



ry rv 



a V(25 kpc) c 



(23) 



with a = 0.17. We studied both annihilating (q oc p 2 ) and decaying (q oc p) dark matter 
scenarios. 



5.1 Flux Ratios 

Figure 2 shows the flux calculated with the exponential diffusion coefficient model described 
above, with V e = 5 km/s, L = 40 kpc, R = 20 kpc, and other parameters as specified in 
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Sec. 4, as a function of antiproton energy. The flux is normalized to the flux calculated 
within conventional formalism (MED propagation model). Results are shown for Einasto 
and isothermal profiles, and for annihilating and decaying dark matter. In all cases, we show 
numerical results along with the analytic approximations (for cases 1 and 2 from Section 4). 
We plot the ratios (rather than actual fluxes) since they are independent of the antiproton 
injection spectrum, dark matter annihilation cross section or decay width, and other such 
parameters, and highlight the variation from adopting the different propagation models. 

Since the diffusion coefficient increases faster with energy in the conventional formalism, 
antiproton flux decreases more rapidly, leading to a steady rise in the ratio of fluxes at high 
energies in all panels of Figure 2. The downward sloping behavior of the ratios at energies 
lower than ~30 GeV is due to the fact that the convective wind term becomes important at 
these energies, where it has a comparatively stronger effect in the conventional formalism. 
In all cases, the calculated antiproton fluxes between 15 to 5000 GeV in the two models are 
quite close to each other, differing at most by about 25%. The difference is smaller than 
the systematic uncertainty, which can be estimated as the difference between the results 
of applying the conventional formalism with the MIN, MED and MAX parameters (see, for 
example, fig. 7 in Ref. [6]). Thus, effects of position-dependent diffusion coefficient do not yet 
need to be included in deriving constraints on dark matter models from the antiproton flux 
measurements as long as the present uncertainties on the galactic propagation parameters 
are fully taken into account. On the other hand, the model with exponentially changing 
diffusion coefficient almost certainly better reflects the correct physics than the conventional 
one, and the solutions we obtained are almost as simple as the conventional ones, especially 
at high energies where the wind term can be ignored and analytic approximations can be 
used. We therefore advocate using it as a "benchmark" propagation model in dark matter 
studies. 



The PAMELA experiment recently measured the cosmic ray antiproton flux and antiproton 
to proton ratio up to ~ 180 GeV [2], extending its earlier results [1]. Since no excess of 
antiprotons over what is expected from conventional astrophysics is evident, these measure- 
ments can be used to place bounds on properties of dark matter; this has been done in 
several earlier works (see, for instance, [6, 13, 14, 15, 16]). As an illustration of our results, 
we recalculate these bounds using the propagation model of Section 4. We focus on a stable 
WIMP dark matter pair-annihilating into W + W~ pairs, such as a wino-like neutralino in 
the MSSM, or the lightest T-odd particle of the Littlest Higgs model with T-parity [17]. 
Antiprotons are produced in hadronic decays of the W bosons. 

We use the solution in Eqs. (20)-(21), and set the convective wind term to zero. The 
source term q(E,x) can be written as 



5.2 Effect on PAMELA Constraints 




(24) 
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where T = E — m p is the kinetic energy of the antiproton, and the fragmentation function 
dNp/dT can be parameterized as [8] 

dN- 

-^ = {pix^+ P2 \\og w xrr i ) (25) 

dx 

where x = T/m x . This parameterization fits the results of the PYTHIA Monte-Carlo code 
[18]. The parameters Pi depend on the WIMP mass m x and, for annihilation into W + W~ , 
can be written as [8] 

pi = (306.0m;- 28 + 7.2 x lO^mJ 25 )" 1 , 
p 2 = (2.32m; 05 )- 1 , 
p 3 = (-8.5m- - 31 )" 1 , 

p 4 = (-0.39m-°- 17 - 2.0 x Kr 2 mJ 23 ) -1 . (26) 

This parameterization is valid for the WIMP mass in the range 50 GeV to 5 TeV. 

To obtain bounds, we fit to the PAMELA antiproton flux data presented in [2] , conserva- 
tively assuming zero antiproton background from astrophysical sources. Figure 3 shows the 
bound on the annihilation cross section (crv), obtained by requiring that the signal not exceed 
the measured flux in any of the energy bins by more than 2a (statistical error only). The 
corresponding constraints calculated within the conventional propagation formalism (MED 
propagation model) are also shown; these are in agreement with similar limits calculated 
in [16]. The model with the exponentially increasing diffusion coefficient is found to relax 
the constraints by ~ 20% for both isothermal and Einasto profiles. In most cases, the bound 
is enforced by the flux in the highest energy bin, at 100-180 GeV; Figure 2 indeed shows 
that there is a difference of about 20% in the antiproton flux calculated in these two models. 

We also plot, for comparison, the theoretically predicted s-wave cross section for the 
MSSM wino annihilating into W + W~ [19] 

. . 1 0.65 4 (l-x) 3 / 2 

M ^m|^V^F' ^ 

where x = m^/m 2 ,. The conventional propagation model rules out wino-like neutralinos 
lighter than ~ 260 GeV (for an isothermal profile); the model with the exponentially in- 
creasing diffusion coefficient slightly weakens this lower bound to ~ 240 GeV. As mentioned 
earlier, these results should be taken with a grain of salt, since in the conventional formalism 
calculation we fixed the parameters corresponding to the MED propagation model; including 
the uncertainties on these parameters would result in a band of predictions encompassing 
those of the model with exponential diffusion coefficient. 

6 Discussion 

In this paper, we extended the conventional approach to calculation of the antiproton flux 
from dark matter annihilation/decay to allow for the possibility of a position-dependent 
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dark matter mass (GeV) dark matter mass (GeV) 



Figure 3: Upper limits on annihilation cross section as a function of m x for XX ~ >* W + W~ for 
isothermal (left) and Einasto (right) profiles. Solid (dot-dashed) curves denote limits from 
using the propagation model with exponentially increasing (constant) diffusion coefficient. 
The theoretical prediction for s-wave LSP wino annihilation is also given (dashed curve). 

diffusion coefficient. We studied two models, the three- zone model where the diffusion co- 
efficient jumps on the boundary between the "diffusion" and "free-propagation" zones, and 
a model with the diffusion coefficient growing smoothly (exponentially) with distance away 
from the disk. In the first model, we found an analytic solution and showed that in the limit 
of infinite diffusion coefficient in the free-propagation zone the flux on Earth is not modified 
by the sources in that zone due to perfect reflection of antiprotons from the zone boundary. 
This seems to justify the conventional calculation even in situations when there are sources 
outside the diffusion zone, as is common in dark matter studies. On the other hand, the 
three-zone model is a rather crude approximation, since magnetic fields, and with them the 
diffusion coefficient, are expected to vary smoothly with distance from the disk. Such smooth 
variation was incorporated in our second model, which assumed an exponentially growing 
diffusion coefficient, which has been shown to produce consistent fits to conventional astro- 
physical cosmic ray observables in Ref. [4]. We found numerical solutions as well as analytic 
approximations valid at high antiproton energies for this model. We found that the resulting 
fluxes on Earth differ from the predictions of the conventional calculation with the MED 
propagation model by at most about 25%. Such deviations are well within the numerous un- 
certainties inherent in the calculation, hence the use of the conventional approach is justified, 
at least at present, in computing the antiproton fluxes from dark matter annihilation/decay 
and using them to place bounds on dark matter models. Nevertheless, since the model with 
exponentially growing diffusion coefficient almost certainly captures the correct physics of 
charged particle propagation in the galaxy better than the conventional one, and since the 
solutions we obtained in this model - especially the analytic solutions, which are accurate 
at high energies - are essentially as simple as the conventional ones, we advocate the use of 
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this model as a benchmark for dark matter studies. 
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